data <- readRDS("very_low_birthweight.RDS")
data <- data %>%
mutate(id = row_number()) %>%
select(id, everything())
data_cleaned <- data %>%
select(where(~ sum(is.na(.)) <= 100)) %>%
drop_na() %>%
mutate(across(c(twn, vent, pneumo, pda, cld, dead), as.factor)) %>%
mutate(hospstay = as.numeric(hospstay))
theme_custom <- theme(
panel.background = element_rect(fill = "white", color = "black"),
panel.grid.major = element_line(color = "grey"),
panel.grid.minor = element_blank(),
axis.title = element_text(size = 15),
axis.text = element_text(size = 12),
axis.line = element_line(color = "black"),
plot.title = element_text(size = 15, face = "bold", hjust = 0.5, color = "#2c2c2c"),
axis.text.x = element_text(size = 14, angle = 60, hjust = 1),
axis.text.y = element_text(size = 14),
strip.text = element_text(size = 16, face = "bold"),
legend.background = element_rect(fill = "white", color = "black"),
legend.key = element_rect(fill = "white"),
legend.position = "top"
)
data_cleaned <- data_cleaned %>%
mutate(across(where(is.numeric), ~ ifelse(
. < quantile(., 0.25, na.rm = TRUE) - 1.5 * IQR(., na.rm = TRUE) |
. > quantile(., 0.75, na.rm = TRUE) + 1.5 * IQR(., na.rm = TRUE),
NA, .
))) %>%
drop_na() %>%
mutate(across(where(is.character), as.factor))
numeric_vars <- data_cleaned %>%
select(where(is.numeric), -id)
numeric_vars %>%
select(where(is.numeric)) %>%
reshape2::melt() %>%
ggplot(aes(x = value, fill = variable)) +
geom_density(alpha = 0.5, colour = "grey49") +
facet_wrap(~ variable, scales = "free") +
labs(
title = "Плотность числовых переменных",
x = "Значение",
y = "Плотность",
fill = "Переменная"
) +
theme_custom
## No id variables; using all as measure variables
data_cleaned %>%
select(inout, bwt, gest) %>%
reshape2::melt(id.vars = "inout") %>%
ggplot(aes(x = value, fill = inout)) +
geom_density(alpha = 0.5) +
facet_wrap(~ variable, scales = "free") +
labs(
title = "Плотности веса при рождении и гестационного возраста по In/Out",
x = "Значение",
y = "Плотность",
fill = "In/Out"
) +
theme_custom
test_results <- data_cleaned %>%
rstatix::t_test(lowph ~ inout, var.equal = FALSE) %>%
add_significance()
ggbarplot(
data_cleaned,
x = "inout",
y = "lowph",
add = c("mean_se", "jitter"),
fill = "inout",
color = "black",
palette = "Dark2",
width = 0.3
) +
geom_hline(
yintercept = mean(data_cleaned$lowph, na.rm = TRUE),
linetype = "dashed",
color = "gray50",
size = 1
) +
stat_pvalue_manual(
test_results,
label = "p = {p}",
y.position = max(data_cleaned$lowph, na.rm = TRUE) + 0.5,
size = 5,
color = "black"
) +
scale_y_continuous(
limits = c(0, max(data_cleaned$lowph, na.rm = TRUE) + 1),
expand = expansion(mult = c(0, 0.05))
) +
labs(
title = "Средние значения lowph по группам",
x = "Группа",
y = "lowph"
) +
theme_custom
Можно сделать вывод, что дети, доставленные в госпиталь из других мест, могут иметь более высокий уровень смертности.
cont_data <- data_cleaned %>%
select(where(is.numeric)) %>%
select(-birth, -year, -exit, -id)
correlation_matrix <- cor(cont_data, use = "complete.obs", method = "spearman")
print(correlation_matrix)
## hospstay lowph pltct bwt gest apg1
## hospstay 1.0000000 -0.2690291 -0.14567930 -0.4781800 -0.41533541 -0.1169665
## lowph -0.2690291 1.0000000 0.28421133 0.3130708 0.35814260 0.2411210
## pltct -0.1456793 0.2842113 1.00000000 0.2122459 0.04150965 0.2833179
## bwt -0.4781800 0.3130708 0.21224590 1.0000000 0.68555612 0.3091558
## gest -0.4153354 0.3581426 0.04150965 0.6855561 1.00000000 0.2466861
## apg1 -0.1169665 0.2411210 0.28331790 0.3091558 0.24668614 1.0000000
pheatmap(
correlation_matrix,
display_numbers = TRUE,
color = colorRampPalette(c("lightblue", "white", "red"))(50),
main = "Тепловая карта корреляционной матрицы"
)
correlation_matrix %>%
network_plot(min_cor = .0,
legend ="full",
colors = c("red", "white", "blue")) +
labs(title = "Сетевой график корреляций") +
theme_custom
scaled_data <- cont_data %>%
scale()
distance_matrix <- dist(scaled_data, method = "euclidean")
hclust_model <- hclust(distance_matrix, method = "ward.D2")
k <- 4
fviz_dend(
hclust_model,
k = k,
rect = TRUE,
rect_border = "jco",
rect_fill = TRUE,
cex = 0.8,
main = "Дендрограмма иерархической кластеризации"
)
pheatmap(
scaled_data,
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
clustering_method = "ward.D2",
color = colorRampPalette(c("blue", "white", "red"))(50),
main = "Heatmap с иерархической кластеризацией",
scale = "row",
display_numbers = FALSE,
fontsize_row = 10,
fontsize_col = 10
)
Интерпретация: Более высокий гестационный возраст (gest) связан с более высоким весом при рождении (bwt), что соответствует медицинским фактам. Недоношенные дети имеют меньший вес, что подтверждается на графике. Длительность госпитализации (hospstay) зависит от комплексного набора факторов, таких как вес при рождении (bwt), гестационный возраст (gest), и состояние ребенка (lowph, apg1). Дети с низким bwt и gest, как правило, госпитализируются на более длительный срок, что подтверждается красными зонами тепловой карты.
pca_data <- prcomp(cont_data,scale = TRUE)
summary(pca_data)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 1.5954 1.0661 0.8823 0.8512 0.73334 0.52677
## Proportion of Variance 0.4242 0.1894 0.1298 0.1208 0.08963 0.04625
## Cumulative Proportion 0.4242 0.6136 0.7434 0.8641 0.95375 1.00000
fviz_eig(pca_data, addlabels = TRUE, ylim = c(0, 50), main = "Доля объяснённой дисперсии")
fviz_contrib(pca_data, choice = "var", axes = 1, top = 24)
fviz_contrib(pca_data, choice = "var", axes = 2, top = 24)
fviz_contrib(pca_data, choice = "var", axes = 3, top = 24)
fviz_pca_var(
pca_data,
col.var = "contrib",
repel = TRUE,
title = "График переменных PCA"
)
fviz_pca_ind(
pca_data,
geom = "point",
col.ind = "cos2",
repel = TRUE,
title = "График наблюдений PCA"
)
Интерпретация: Первая компонента (Dim1) объясняет 41.8% дисперсии данных, вторая — 18.9%, третья — 13.2%. В сумме первые три компоненты объясняют 73.9% общей дисперсии, что свидетельствует о том, что они хорошо представляют исходные данные. Значительное снижение дисперсии после третьей компоненты говорит о том, что для анализа можно сосредоточиться на первых трёх компонентах.
Переменные bwt (вес при рождении) и gest (гестационный возраст) вносят наибольший вклад в первую компоненту. Это указывает на то, что первая компонента в основном отражает характеристики, связанные с физическими параметрами при рождении.
Вторая компонента в значительной степени определяется переменной pltct (количество тромбоцитов), которая вносит около 50% вклада. Она отражает характеристики, связанные с состоянием здоровья.
Третья компонента определяется переменной apg1 (оценка по шкале Апгар на первой минуте), которая вносит около 50% вклада. Это может свидетельствовать о том, что она отражает показатели здоровья ребенка при рождении.
Судя по графику переменных PCA bwt и gest имеют близкие направления, что подтверждает их положительную корреляцию. pltct и apg1 имеют небольшую корреляцию с другими переменными.
Применение шкалирования было необходимо, так как переменные имели разные единицы измерения (например, вес в граммах и тромбоциты в единицах).
ggbiplot(pca_data,
scale=0,
groups = data_cleaned$dead,
ellipse = T,
alpha = 0.4) +
scale_colour_brewer(palette = "Set2") +
labs(
title = "PCA Biplot",
) +
theme_custom
data_id <- cont_data %>%
mutate(id = as.integer(row_number()))
biplot <- ggbiplot(
pca_data,
scale = 0,
groups = data_cleaned$dead,
ellipse = TRUE,
alpha = 0.4
) +
geom_point(
aes(text = paste("ID:", data_id$id)),
alpha = 0
) +
scale_colour_brewer(palette = "Set2") +
labs(
title = "PCA Biplot: Интерактивный анализ",
color = "Статус 'Dead'"
) +
theme_custom
interactive_biplot <- ggplotly(
biplot,
tooltip = "text"
)
n_cases <- length(unique(data_cleaned$dead))
for (i in 1:n_cases) {
interactive_biplot$x$data[[i]]$name <- ifelse(i == 1, "Alive", "Dead")
interactive_biplot$x$data[[i]]$legendgroup <- i
interactive_biplot$x$data[[i + n_cases]]$name <- ifelse(i == 1, "Alive", "Dead")
interactive_biplot$x$data[[i + n_cases]]$legendgroup <- i
interactive_biplot$x$data[[i + n_cases]]$showlegend <- FALSE
}
interactive_biplot$x$layout$legend$title$text <- "Статус"
interactive_biplot
Интерпретация: PCA-анализ показал, что первые две компоненты объясняют 60.7% вариации данных. PC1 отражает физические параметры, такие как вес при рождении (bwt) и гестационный возраст (gest), а PC2 — длительность госпитализации (hospstay). Эллипсы групп dead сильно перекрываются, что указывает на слабую способность PCA разделять выживших и умерших. Использование колонки dead для выводов некорректно, так как PCA не используется для выявления причинно-следственных связей или ассоциаций между переменными.
umap_config <- umap.defaults
umap_config$n_neighbors <- 15
umap_config$min_dist <- 0.1
umap_config$metric <- "euclidean"
umap_result <- umap(scaled_data, config = umap_config)
umap_data <- as.data.frame(umap_result$layout)
colnames(umap_data) <- c("UMAP1", "UMAP2")
umap_data$dead <- data_cleaned$dead
umap_plot <- ggplot(umap_data, aes(x = UMAP1, y = UMAP2, color = dead)) +
geom_point(alpha = 0.7, size = 3) +
scale_color_brewer(palette = "Set2") +
labs(
title = "Проекция данных в 2D-пространство",
x = "UMAP1",
y = "UMAP2",
color = "Статус 'Dead'"
) +
theme_custom
umap_plot
Интерпретация: PCA показывает глобальную структуру данных с учетом дисперсии, тогда как UMAP подчеркивает локальные зависимости. UMAP удобен для анализа сходства между наблюдениями, но не предназначен для изучения связей между переменными
umap_variations <- function(data, n_neighbors, min_dist, target_column) {
umap_data <- recipe(~., data = data) %>%
step_normalize(all_predictors()) %>%
step_umap(all_predictors(), num_comp = 2, neighbors = n_neighbors, min_dist = min_dist) %>%
prep() %>%
juice()
umap_data <- umap_data %>%
mutate(dead = target_column)
ggplot(umap_data, aes(UMAP1, UMAP2)) +
geom_point(aes(color = as.factor(dead)), alpha = 0.6, size = 2) +
scale_color_brewer(palette = "Set2") +
labs(
title = paste("n_neighbors =", n_neighbors, ", min_dist =", min_dist),
color = "Dead"
) +
theme_custom
}
plot1 <- umap_variations(cont_data, n_neighbors = 10, min_dist = 0.01, target_column = data_cleaned$dead)
plot2 <- umap_variations(cont_data, n_neighbors = 50, min_dist = 0.01, target_column = data_cleaned$dead)
plot3 <- umap_variations(cont_data, n_neighbors = 10, min_dist = 0.1, target_column = data_cleaned$dead)
plot4 <- umap_variations(cont_data, n_neighbors = 50, min_dist = 0.1, target_column = data_cleaned$dead)
(plot1 | plot2) / (plot3 | plot4)
Интерпретация: При меньших значениях n_neighbors, алгоритм больше сосредотачивается на локальных связях между точками, что приводит к более разрозненным кластерам. При больших же значениях, алгоритм учитывает более глобальную структуру, сглаживая распределение данных.
Меньшие значения min_dist способствуют более плотным кластерам, что подчеркивает локальную структуру данных. Более высокие значения приводят к более равномерному распределению точек, что снижает визуальную плотность кластеров.
data_50perm <- cont_data %>%
mutate(bwt = case_when(
row_number() %in% sample(row_number(), size = n() * 0.5) ~ sample(bwt),
TRUE ~ bwt
))
data_100perm <- cont_data %>%
mutate(bwt = sample(bwt))
cont_data_50perm <- data_50perm %>%
select(where(is.numeric))
cont_data_100perm <- data_100perm %>%
select(where(is.numeric))
pca_50perm <- prcomp(cont_data_50perm, scale = TRUE)
pca_100perm <- prcomp(cont_data_100perm, scale = TRUE)
umap_50perm <- umap(scale(cont_data_50perm))
umap_100perm <- umap(scale(cont_data_100perm))
umap_data_50 <- as.data.frame(umap_50perm$layout)
colnames(umap_data_50) <- c("UMAP1", "UMAP2")
umap_data_50$dead <- data_cleaned$dead
plot_umap_50 <- ggplot(umap_data_50, aes(x = UMAP1, y = UMAP2, color = dead)) +
geom_point(alpha = 0.7, size = 3) +
scale_color_brewer(palette = "Set2") +
labs(
title = "50% Permuted Data",
x = "UMAP1",
y = "UMAP2",
color = "Dead"
) +
theme_custom
umap_data_100 <- as.data.frame(umap_100perm$layout)
colnames(umap_data_100) <- c("UMAP1", "UMAP2")
umap_data_100$dead <- data_cleaned$dead
plot_umap_100 <- ggplot(umap_data_100, aes(x = UMAP1, y = UMAP2, color = dead)) +
geom_point(alpha = 0.7, size = 3) +
scale_color_brewer(palette = "Set2") +
labs(
title = "100% Permuted Data",
x = "UMAP1",
y = "UMAP2",
color = "Dead"
) +
theme_custom
(plot_umap_50 | plot_umap_100)
biplot_50 <- ggbiplot(pca_50perm,
scale = 0,
alpha = 0.4,
ellipse = TRUE,
groups = data_cleaned$dead) +
labs(
title = "PCA (50% permutated)",
fill = "dead",
color = "dead"
) +
theme_custom
biplot_100 <- ggbiplot(pca_100perm,
scale = 0,
alpha = 0.4,
ellipse = TRUE,
groups = data_cleaned$dead) +
labs(
title = "PCA (100% permutated)",
fill = "dead",
color = "dead"
) +
theme_custom
(biplot_50 | biplot_100)
Интерпретация: Пермутация данных значительно влияет на результаты как
PCA, так и UMAP. В PCA кумулятивный процент объяснённой вариации заметно
уменьшается: при 50% пермутации структура данных частично сохраняется,
но кластеризация ослаблена, а при 100% данные становятся случайными, и
кластеры исчезают. UMAP, напротив, более устойчив к частичной
пермутации, сохраняя локальные связи при 50% изменении данных, но при
100% пермутации распределение точек становится хаотичным, аналогично
PCA. Визуализация UMAP лучше сохраняет локальную структуру при частичной
пермутации, в то время как PCA демонстрирует потерю глобальной структуры
уже на этапе 50% изменений.
data_cleaned_imp <- data %>%
select(hospstay, lowph, pltct, bwt, gest, apg1, dead) %>%
mutate(across(
where(is.numeric),
~ ifelse(is.na(.), mean(., na.rm = TRUE), .)
)) %>%
mutate(
hospstay = as.numeric(hospstay),
dead = as.factor(dead))
numeric_imp <- data_cleaned_imp %>%
select(hospstay, lowph, pltct, bwt, gest, apg1)
correlation_matrix_imput <- cor(numeric_imp, use = "complete.obs", method = "spearman")
pheatmap(
correlation_matrix_imput,
display_numbers = TRUE,
color = colorRampPalette(c("lightblue", "white", "red"))(50),
main = "Тепловая карта корреляционной матрицы imp"
)
correlation_matrix_imput %>%
network_plot(min_cor = .0,
legend ="full",
colors = c("red", "white", "blue")) +
labs(title = "Сетевой график корреляций imp") +
theme_custom
scaled_data_imput <- numeric_imp %>%
scale() %>%
as.data.frame() %>%
filter(apply(., 1, var) > 0)
distance_matrix <- dist(scaled_data_imput, method = "euclidean")
hclust_model <- hclust(distance_matrix, method = "ward.D2")
k <- 4
fviz_dend(
hclust_model,
k = k,
rect = TRUE,
rect_border = "jco",
rect_fill = TRUE,
cex = 0.8,
main = "Дендрограмма иерархической кластеризации imp"
)
pheatmap(
scaled_data_imput,
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
clustering_method = "ward.D2",
color = colorRampPalette(c("blue", "white", "red"))(50),
main = "Heatmap с иерархической кластеризацией imp",
scale = "row",
display_numbers = FALSE,
fontsize_row = 10,
fontsize_col = 10
)
Интерпретация: графики с импутацией выглядят более равномерными, тогда
как графики оригинальных данных лучше отражают реальные зависимости.
Оригинальный подход с удалением пропусков обеспечивает более точные и
достоверные результаты, так как анализируется только полная информация.
Однако он теряет часть данных, что может снижать репрезентативность при
большом количестве пропусков. Импутация, напротив, сохраняет полный
объем данных, что делает визуализации более полными, но может сглаживать
вариативность и создавать ложные взаимосвязи.
data_for_pca_imp <- data_cleaned_imp %>%
select(-dead) %>%
scale()
pca_data_imp <- prcomp(data_for_pca_imp, scale = TRUE)
fviz_eig(pca_data_imp, addlabels = TRUE, ylim = c(0, 50), main = "Доля объяснённой дисперсии imp")
fviz_contrib(pca_data_imp, choice = "var", axes = 1, top = 24)
fviz_contrib(pca_data_imp, choice = "var", axes = 2, top = 24)
fviz_contrib(pca_data_imp, choice = "var", axes = 3, top = 24)
fviz_pca_var(
pca_data_imp,
col.var = "contrib",
repel = TRUE,
title = "График переменных PCA imp"
)
fviz_pca_ind(
pca_data_imp,
geom = "point",
col.ind = "cos2",
repel = TRUE,
title = "График наблюдений PCA imp"
)
ggbiplot(pca_data_imp,
scale=0,
groups = data_cleaned_imp$dead,
ellipse = T,
alpha = 0.4) +
scale_colour_brewer(palette = "Set2") +
labs(
title = "PCA Biplot imp",
) +
theme_custom
umap_result_imp <- umap(data_for_pca_imp, config = umap_config)
umap_data_imp <- as.data.frame(umap_result_imp$layout)
colnames(umap_data_imp) <- c("UMAP1", "UMAP2")
umap_data_imp$dead <- data_cleaned_imp$dead
ggplot(umap_data_imp, aes(x = UMAP1, y = UMAP2, color = dead)) +
geom_point(alpha = 0.7, size = 3) +
scale_color_brewer(palette = "Set2") +
labs(
title = "Проекция данных в 2D-пространство imp",
x = "UMAP1",
y = "UMAP2",
color = "Статус 'Dead'"
) +
theme_custom
Интерпретация: Графики с импутацией демонстрируют смещение вклада переменных, где hospstay доминирует в третьей главной компоненте, объясняя более 95% ее вклада. Это приводит к изменению общей структуры данных, делая кластеры с умершими детьми более явными, особенно на PCA-биплоте. Это также проявляется в меньшем проценте объясненной дисперсии для первой компоненты (39.4% против 42.4% в оригинальных данных) и второй компоненты (17.2% против 18.9%), что указывает на более равномерное распределение дисперсии между компонентами при импутации. Визуализация UMAP показывает лучшее разделение групп, однако сглаживание из-за импутации делает общую структуру менее точной.